2.1 Basic Concepts

The basic ideas of association rule mining and basket analysis draws on basic ideas from probability theory.

We speak in terms of the itemset: that is, a collection of one or more items that co-occur in a transaction.

For example, suppose we have a list of transactions as follows:

Using the above set of transactions, and itemset may be “milk” or “bread, butter”.

The support of an itemset \(X\), \(\text{Supp}(X)\), is defined as the proportion of transactions in the dataset which contain the itemset.

In the above example:

\[ \text{Supp}(\text{\{milk, bread\}}) = \frac{2}{5} = 0.40. \]

A rule, \(X \Rightarrow Y\), between two itemsets \(X\) and \(Y\) is a directed relationship of the itemset \(X\) showing the presence of \(Y\). The rule is not symmetric: \(X \Rightarrow Y\) and \(Y \Rightarrow X\) are not the same.

The confidence for the rule \(X \implies Y\), \(\text{Conf}(X \Rightarrow Y)\) is defined by

\[ \text{Conf}(X \Rightarrow Y) = \frac{\text{Supp}(X \cup Y)}{\text{Supp}(X)}. \]

So, to calculate the confidence for a rule:

\[ \text{Conf}(\text{\{milk, bread\}} \Rightarrow \text{\{butter\}}) = \frac{0.2}{0.4} = 0.5. \]

To illustrate how rules are not symmetric:

\[ \text{Conf}(\text{\{butter\}} \Rightarrow \text{\{milk, bread\}}) = \frac{0.2}{0.6} = 0.33. \]

Finally, we want a measure of the strength of the relationship between the itemsets \(X\) and \(Y\). That is, measuring the effect of the presence of \(X\) on the presence of \(Y\). We measure this by defining the lift of a rule as

\[ \text{Lift}(X \Rightarrow Y) = \frac{\text{Supp}(X \cup Y)}{\text{Supp}(X) \text{Supp}(Y)}. \]

Again, we repeat our calculations for our rule.

\[ \text{Lift}(\text{\{bread, milk\}} \Rightarrow \text{\{butter\}}) = \frac{0.2}{(0.4)(0.6)} = \frac{0.2}{0.24} = 0.8333 \]

Lift values greater than 1 implies the presence of \(X\) increases the probability of \(Y\) being present when compared to the unconditional probability.

Now that we have these metrics and concepts, we can turn our attention to trying to find rules in a given dataset, using these metrics to rank them.

Rather than using brute-force approaches to discovering these rules, we use a number of different algorithms to find associations within the dataset.

The two main algorithms for discovering some rules are the apriori and the eclat algorithms.

2.2 Construct apriori Rules

We now want to construct the association rules using the apriori algorithm. To do this, we need to set parameters such as the minimum support and the minimum confidence level.

This gives us a set of association rules, along with the support and lift.

basket_apriori <- apriori(
    basket_arules,
    parameter = list(supp = 0.005, conf = 0.8)
    )

## Apriori
## 
## Parameter specification:
##  confidence minval smax arem  aval originalSupport maxtime support minlen
##         0.8    0.1    1 none FALSE            TRUE       5   0.005      1
##  maxlen target  ext
##      10  rules TRUE
## 
## Algorithmic control:
##  filter tree heap memopt load sort verbose
##     0.1 TRUE TRUE  FALSE TRUE    2    TRUE
## 
## Absolute minimum support count: 197 
## 
## set item appearances ...[0 item(s)] done [0.00s].
## set transactions ...[4895 item(s), 39510 transaction(s)] done [0.34s].
## sorting and recoding items ... [1445 item(s)] done [0.03s].
## creating transaction tree ... done [0.02s].
## checking subsets of size 1 2 3 4 5 6 done [0.19s].
## writing ... [537 rule(s)] done [0.00s].
## creating S4 object  ... done [0.01s].

basket_apriori_tbl <- basket_apriori %>%
  as("data.frame") %>%
  as_tibble() %>%
  arrange(desc(lift))

basket_apriori_tbl %>% glimpse()

## Rows: 537
## Columns: 6
## $ rules      <chr> "{22917,22918} => {22916}", "{22916,22919} => {22917}", "{2…
## $ support    <dbl> 0.005391040, 0.005188560, 0.005391040, 0.005087320, 0.00528…
## $ confidence <dbl> 0.9466667, 0.9579439, 0.9508929, 0.9392523, 0.9500000, 0.93…
## $ coverage   <dbl> 0.005694761, 0.005416350, 0.005669451, 0.005416350, 0.00556…
## $ lift       <dbl> 150.2120, 150.1919, 149.0864, 149.0356, 148.9464, 148.7131,…
## $ count      <int> 213, 205, 213, 201, 209, 209, 205, 204, 204, 203, 203, 232,…

We now want to inspect this table using the ruleExplorer()

basket_apriori %>% ruleExplorer()

To help visualise these rules, we can produce a basic scatterplot of the metrics.

ggplot(basket_apriori_tbl) +
  geom_point(aes(x = confidence, y = lift), alpha = 0.2) +
  xlab("Rule Confidence") +
  ylab("Rule Lift") +
  ggtitle("Scatterplot of Association Rule Metrics")

2.3 Construct eclat Rules

An alternative method of constructing association rules is to use the eclat algorithm. The code for doing this is slightly different, but gives us similar outputs.

basket_eclat <- eclat(
    basket_arules,
    parameter = list(support = 0.005)
    ) %>%
  ruleInduction(
    basket_arules,
    confidence = 0.8
    )

## Eclat
## 
## parameter specification:
##  tidLists support minlen maxlen            target  ext
##     FALSE   0.005      1     10 frequent itemsets TRUE
## 
## algorithmic control:
##  sparse sort verbose
##       7   -2    TRUE
## 
## Absolute minimum support count: 197 
## 
## create itemset ... 
## set transactions ...[4895 item(s), 39510 transaction(s)] done [0.33s].
## sorting and recoding items ... [1445 item(s)] done [0.02s].
## creating sparse bit matrix ... [1445 row(s), 39510 column(s)] done [0.01s].
## writing  ... [7742 set(s)] done [4.12s].
## Creating S4 object  ... done [0.00s].

basket_eclat_tbl <- basket_eclat %>%
  as("data.frame") %>%
  as_tibble() %>%
  arrange(desc(lift))

basket_eclat_tbl %>% glimpse()

## Rows: 537
## Columns: 5
## $ rules      <chr> "{22917,22918} => {22916}", "{22916,22919} => {22917}", "{2…
## $ support    <dbl> 0.005391040, 0.005188560, 0.005391040, 0.005087320, 0.00528…
## $ confidence <dbl> 0.9466667, 0.9579439, 0.9508929, 0.9392523, 0.9500000, 0.93…
## $ lift       <dbl> 150.2120, 150.1919, 149.0864, 149.0356, 148.9464, 148.7131,…
## $ itemset    <int> 69, 41, 69, 39, 73, 73, 41, 26, 26, 38, 68, 76, 76, 40, 69,…

Once again, we inspect the data using ruleExplorer()

basket_eclat %>% ruleExplorer()

2.4 Compare Algorithms

We now want to compare the outputs of both algorithms in terms of association rules and how they compare.

basket_ap_tbl <- basket_apriori_tbl %>%
  select(rules, support, confidence, lift)

basket_ec_tbl <- basket_eclat_tbl %>%
  select(rules, support, confidence, lift)

rules_comparison_tbl <- basket_ap_tbl %>%
  full_join(basket_ec_tbl, by = "rules", suffix = c("_a", "_e"))

rules_comparison_tbl %>% glimpse()

## Rows: 537
## Columns: 7
## $ rules        <chr> "{22917,22918} => {22916}", "{22916,22919} => {22917}", "…
## $ support_a    <dbl> 0.005391040, 0.005188560, 0.005391040, 0.005087320, 0.005…
## $ confidence_a <dbl> 0.9466667, 0.9579439, 0.9508929, 0.9392523, 0.9500000, 0.…
## $ lift_a       <dbl> 150.2120, 150.1919, 149.0864, 149.0356, 148.9464, 148.713…
## $ support_e    <dbl> 0.005391040, 0.005188560, 0.005391040, 0.005087320, 0.005…
## $ confidence_e <dbl> 0.9466667, 0.9579439, 0.9508929, 0.9392523, 0.9500000, 0.…
## $ lift_e       <dbl> 150.2120, 150.1919, 149.0864, 149.0356, 148.9464, 148.713…

2.5 Reducing Minimum Confidence

While high confidence rules are useful, they are more likely to find rules that are “obvious” as the probabilities are such that co-occuring basket items will be noticed as being together - or possibly be natural complements: butter, milk and bread is a good example.

Instead, we are also interested in less obvious rules, and so we reduce our confidence threshold and see how many additional rules are discovered.

basket_lower_rules <- apriori(
    basket_arules,
    parameter = list(supp = 0.005, conf = 0.4)
  )

## Apriori
## 
## Parameter specification:
##  confidence minval smax arem  aval originalSupport maxtime support minlen
##         0.4    0.1    1 none FALSE            TRUE       5   0.005      1
##  maxlen target  ext
##      10  rules TRUE
## 
## Algorithmic control:
##  filter tree heap memopt load sort verbose
##     0.1 TRUE TRUE  FALSE TRUE    2    TRUE
## 
## Absolute minimum support count: 197 
## 
## set item appearances ...[0 item(s)] done [0.00s].
## set transactions ...[4895 item(s), 39510 transaction(s)] done [0.33s].
## sorting and recoding items ... [1445 item(s)] done [0.03s].
## creating transaction tree ... done [0.02s].
## checking subsets of size 1 2 3 4 5 6 done [0.18s].
## writing ... [5950 rule(s)] done [0.01s].
## creating S4 object  ... done [0.01s].

basket_lower_rules_tbl <- basket_lower_rules %>%
  as("data.frame") %>%
  as_tibble() %>%
  arrange(desc(lift))

ggplot(basket_lower_rules_tbl) +
  geom_point(aes(x = confidence, y = lift), alpha = 0.2) +
  xlab("Rule Confidence") +
  ylab("Rule Lift") +
  ggtitle("Scatterplot of Association Rule Metrics")

3.1 Extract Connected Product Labels

First we want to look at the different disjoint components of the graph, and label them with an ID.

apriori_rules_tblgraph <- apriori_rules_igraph %>%
  igraph::as.undirected(mode = "collapse") %>%
  as_tbl_graph() %>%
  mutate(
    component_id = group_components()
    ) %>%
  group_by(component_id) %>%
  mutate(
    component_size = n()
    ) %>%
  ungroup()

We then want to create groups of common products that form a disjoint cluster within this graph.

product_groups_all_tbl <- apriori_rules_tblgraph %>%
  activate(nodes) %>%
  as_tibble() %>%
  filter(are_na(support)) %>%
  group_by(component_id) %>%
  mutate(
    product_count = n()
    ) %>%
  ungroup() %>%
  select(
    product_group_id = component_id,
    product_count,
    stock_code = label
    ) %>%
  arrange(product_group_id, stock_code)

product_groups_all_tbl %>% glimpse()

## Rows: 142
## Columns: 3
## $ product_group_id <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
## $ product_count    <int> 67, 67, 67, 67, 67, 67, 67, 67, 67, 67, 67, 67, 67, 6…
## $ stock_code       <chr> "20711", "20712", "20713", "20718", "20719", "20723",…

For display purposes, we can show all the stock_id values in a list.

apriori_rules_large_tblgraph <- apriori_rules_tblgraph %>%
  to_subgraph(component_size == max(component_size)) %>%
  use_series(subgraph) %>%
  morph(to_undirected) %>%
  mutate(
    sub_id = group_louvain()
    ) %>%
  unmorph()

product_groups_largest_tbl <- apriori_rules_large_tblgraph %>%
  activate(nodes) %>%
  as_tibble() %>%
  filter(are_na(support)) %>%
  group_by(sub_id) %>%
  mutate(
    product_count = n()
    ) %>%
  ungroup() %>%
  select(
    product_group_id = sub_id, product_count, stock_code = label
    ) %>%
  arrange(product_group_id, stock_code)

product_groups_largest_tbl %>% glimpse()

## Rows: 67
## Columns: 3
## $ product_group_id <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2,…
## $ product_count    <int> 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 16, 1…
## $ stock_code       <chr> "20719", "20723", "20724", "20725", "20727", "20728",…

apriori_rules_large_tblgraph %>%
  toVisNetworkData(idToLabel = FALSE) %>%
  visNetwork(
    nodes = .$nodes %>% transmute(id, label, group = sub_id),
    edges = .$edges
    )

3.2 Evaluating Product Groups

How do we go about assessing the validity of these product groups?

Note that this work is exploratory - in effect this is more sophisticated data exploration. Rather than use this model to make predictions - a job we will need to do at some point, we instead just want to assess how novel these grouping are.

To that end, it may be useful to check the co-occurrence of these products as a group - for each purchase made by a customer, what proportion of the group was featured in this data?

This question is worth exploring, so we should write some code to assess this.

Before we do this, we combine our two lists of product groups into a single table.

product_groups_tbl <- list(
    ALL = product_groups_all_tbl,
    LRG = product_groups_largest_tbl
    ) %>%
  bind_rows(.id = "type") %>%
  mutate(
    group_label = sprintf("%s_%02d", type, product_group_id)
    ) %>%
  group_by(group_label) %>%
  mutate(
    group_size  = n()
    ) %>%
  ungroup() %>%
  select(group_label, group_size, stock_code)

product_groups_tbl %>% glimpse()

## Rows: 209
## Columns: 3
## $ group_label <chr> "ALL_01", "ALL_01", "ALL_01", "ALL_01", "ALL_01", "ALL_01"…
## $ group_size  <int> 67, 67, 67, 67, 67, 67, 67, 67, 67, 67, 67, 67, 67, 67, 67…
## $ stock_code  <chr> "20711", "20712", "20713", "20718", "20719", "20723", "207…

tnx_groups_tbl <- tnx_data_tbl %>%
  select(invoice_id, invoice_date, stock_code) %>%
  group_nest(invoice_id, .key = "invoice_data")

group_props_tbl <- product_groups_tbl %>%
  group_nest(group_label, group_size, .key = "stock_data") %>%
  filter(group_size > 1, group_size < 15) %>%
  expand_grid(tnx_groups_tbl) %>%
  mutate(
    comb_data = future_map2(
      invoice_data, stock_data,
      inner_join,
      by = "stock_code",
    
      .options = furrr_options(globals = FALSE)
      ),
    match_count = map_int(comb_data, nrow),
    group_prop  = match_count / group_size
    ) %>%
  select(group_label, group_size, group_prop) %>%
  filter(group_prop > 0)

group_props_tbl %>% glimpse()

## Rows: 54,147
## Columns: 3
## $ group_label <chr> "ALL_02", "ALL_02", "ALL_02", "ALL_02", "ALL_02", "ALL_02"…
## $ group_size  <int> 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10…
## $ group_prop  <dbl> 0.2, 0.6, 0.3, 0.3, 0.1, 0.7, 0.2, 0.4, 0.2, 0.4, 0.1, 0.2…

We now create a histogram of the proportions for each group, and this gives us a gauge of the ‘novelty’ of each of these groups.

plot_tbl <- group_props_tbl %>%
  mutate(label = glue("{group_label} ({group_size})"))

ggplot(plot_tbl) +
  geom_histogram(aes(x = group_prop), binwidth = 0.1) +
  facet_wrap(vars(label), scales = "free_y") +
  scale_y_continuous(labels = label_comma()) +
  xlab("Proportion") +
  ylab("Purchase Count") +
  ggtitle("Facetted Histograms of Group Coverages by Product Grouping") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5))

product_groups_tbl %>%
  filter(group_size > 1, group_size < 15) %>%
  mutate(stock_code = stock_code %>% str_trim() %>% str_to_upper()) %>%
  left_join(stock_description_tbl, by = "stock_code") %>%
  datatable()

4.1 Determine Distinct Rules Subgraphs

Having converted the association rules to the graph, we then look at the distinct components of this graph and use these as our first pass at these clusters.

apriori_lower_rules_tblgraph <- apriori_lower_rules_igraph %>%
  igraph::as.undirected(mode = "collapse") %>%
  as_tbl_graph() %>%
  mutate(
    component_id = group_components()
    ) %>%
  group_by(component_id) %>%
  mutate(
    component_size = n()
    ) %>%
  ungroup()

apriori_lower_rules_tblgraph %>% print()

## # A tbl_graph: 5460 nodes and 15770 edges
## #
## # A bipartite simple graph with 88 components
## #
## # Node Data: 5,460 × 11 (active)
##   name  label  type support confidence coverage  lift count order component_id
##   <chr> <chr> <dbl>   <dbl>      <dbl>    <dbl> <dbl> <int> <int>        <int>
## 1 25    1505…     1      NA         NA       NA    NA    NA    NA           14
## 2 27    1505…     1      NA         NA       NA    NA    NA    NA           14
## 3 29    1505…     1      NA         NA       NA    NA    NA    NA           14
## 4 217   20674     1      NA         NA       NA    NA    NA    NA            2
## 5 218   20675     1      NA         NA       NA    NA    NA    NA            2
## 6 219   20676     1      NA         NA       NA    NA    NA    NA            2
## # … with 5,454 more rows, and 1 more variable: component_size <int>
## #
## # Edge Data: 15,770 × 2
##    from    to
##   <int> <int>
## 1    11   461
## 2   217   461
## 3   218   461
## # … with 15,767 more rows

product_groups_lower_all_tbl <- apriori_lower_rules_tblgraph %>%
  activate(nodes) %>%
  as_tibble() %>%
  filter(are_na(support)) %>%
  group_by(component_id) %>%
  mutate(
    product_count = n()
    ) %>%
  ungroup() %>%
  select(product_group_id = component_id, product_count, stock_code = label) %>%
  arrange(product_group_id, stock_code)

product_groups_lower_all_tbl %>% glimpse()

## Rows: 460
## Columns: 3
## $ product_group_id <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
## $ product_count    <int> 216, 216, 216, 216, 216, 216, 216, 216, 216, 216, 216…
## $ stock_code       <chr> "20711", "20712", "20713", "20717", "20718", "20719",…

apriori_lower_rules_bigcomp_tblgraph <- apriori_lower_rules_tblgraph %>%
  to_subgraph(component_size == max(component_size)) %>%
  use_series(subgraph) %>%
  mutate(
    sub_id = group_louvain()
    )

apriori_lower_rules_bigcomp_tblgraph %>% print()

## # A tbl_graph: 4456 nodes and 13807 edges
## #
## # A bipartite simple graph with 1 component
## #
## # Node Data: 4,456 × 12 (active)
##   name  label  type support confidence coverage  lift count order component_id
##   <chr> <chr> <dbl>   <dbl>      <dbl>    <dbl> <dbl> <int> <int>        <int>
## 1 247   20711     1      NA         NA       NA    NA    NA    NA            1
## 2 248   20712     1      NA         NA       NA    NA    NA    NA            1
## 3 249   20713     1      NA         NA       NA    NA    NA    NA            1
## 4 253   20717     1      NA         NA       NA    NA    NA    NA            1
## 5 254   20718     1      NA         NA       NA    NA    NA    NA            1
## 6 255   20719     1      NA         NA       NA    NA    NA    NA            1
## # … with 4,450 more rows, and 2 more variables: component_size <int>,
## #   sub_id <int>
## #
## # Edge Data: 13,807 × 2
##    from    to
##   <int> <int>
## 1     3   217
## 2    98   217
## 3    99   217
## # … with 13,804 more rows

product_groups_lower_bigcomp_tbl <- apriori_lower_rules_bigcomp_tblgraph %>%
  activate(nodes) %>%
  as_tibble() %>%
  filter(are_na(support)) %>%
  group_by(sub_id) %>%
  mutate(
    product_count = n()
    ) %>%
  ungroup() %>%
  select(product_group_id = sub_id, product_count, stock_code = label) %>%
  arrange(product_group_id, stock_code)

product_groups_lower_bigcomp_tbl %>% glimpse()

## Rows: 216
## Columns: 3
## $ product_group_id <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
## $ product_count    <int> 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 2…
## $ stock_code       <chr> "20711", "20712", "20713", "20717", "20718", "21033",…

product_groups_lower_tbl <- list(
    ALL = product_groups_lower_all_tbl,
    LRG = product_groups_lower_bigcomp_tbl
    ) %>%
  bind_rows(.id = "type") %>%
  mutate(
    group_label = sprintf("%s_%02d", type, product_group_id)
    ) %>%
  group_by(group_label) %>%
  mutate(
    group_size = n()
    ) %>%
  ungroup() %>%
  select(group_label, group_size, stock_code)

product_groups_lower_tbl %>% glimpse()

## Rows: 676
## Columns: 3
## $ group_label <chr> "ALL_01", "ALL_01", "ALL_01", "ALL_01", "ALL_01", "ALL_01"…
## $ group_size  <int> 216, 216, 216, 216, 216, 216, 216, 216, 216, 216, 216, 216…
## $ stock_code  <chr> "20711", "20712", "20713", "20717", "20718", "20719", "207…

Now that we have our groups we add on the description fields so that interpretation of the different groupings is easier.

product_groups_lower_tbl %>%
  filter(group_size > 1, group_size != max(group_size)) %>%
  mutate(stock_code = stock_code %>% str_trim() %>% str_to_upper()) %>%
  left_join(stock_description_tbl, by = "stock_code") %>%
  datatable()